from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm


plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.rc('font', size=10)
plt.rc('font', family='SimHei')

m_A = 4866
m_B = 2433
f = 4890
delta_m = 1165.992
omiga = 2.2143
k_3 = 167.8395  # 兴波
k_5 = 80000  # 弹簧劲度系数
k_6 = 10000  # 阻尼系数
ro = 1025
g = 9.8
s = np.pi


def P2(k, alpha):
    # 定义一个方程组（微分方程组）
    def pfun(y, t):
        y0, y1, y2, y3 = y
        return np.array([
            1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k * (
                abs(y2 - y0)) ** alpha * (y2 - y0)),
            y0,
            1 / m_B * (-k_5 * (y3 - y1) - k * (abs(y2 - y0)) ** alpha * (y2 - y0)),
            y2
        ])
    t = np.arange(0, 500, 0.2)  # 创建自变量序列
    soli = odeint(pfun, [0, 0, 0, 0], t)
    re = 0
    for i in range(200, len(soli[:, 2])):
        re += k * (abs(soli[:, 0][i] - soli[:, 2][i])) ** alpha * (soli[:, 0][i] - soli[:, 2][i]) ** 2 * 0.2
    return re / 460


# 创建数据点 - 增加分辨率以获得更光滑的表面
alpha_values = np.linspace(0, 100000, 10)  # 从0到100000，500个点（更密集）
beta_values = np.linspace(0, 1, 10)        # 从0到1，200个点

# 生成网格坐标
Alpha, Beta = np.meshgrid(alpha_values, beta_values, indexing='ij')

#计算Z值
Z = []
for i in tqdm(alpha_values):
    for j in beta_values:
        Z.append(P2(i, j))

# 创建图形
plt.figure(figsize=(8, 6))
# 绘制散点图，使用颜色映射
scatter = plt.scatter(Alpha, Beta, c=Z, cmap='viridis', s=2000, alpha=0.8)
# 添加颜色条
plt.colorbar(scatter, label='值')
# 添加标题和轴标签
plt.title('平均输出功率与比例系数、幂指数曲线')
plt.xlabel('比例系数')
plt.ylabel('幂指数')
# 显示图形
plt.show()
